package scu.maqiang.fes;

import scu.maqiang.mesh.MatrixFunc;
import scu.maqiang.mesh.Mesh3H8;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

import java.util.Arrays;

import static scu.maqiang.numeric.Constants.GP2;
import static scu.maqiang.numeric.Constants.GW2;
import static scu.maqiang.numeric.MVO.*;
import static scu.maqiang.numeric.MVO.matmul;

public class FES3H8PE extends FES{

	public FES3H8PE(Mesh3H8 mesh) {
		super(mesh, 4);
	}

	/**
	 * 计算单元压电刚阵, 注意排列顺序为11, 22, 33, 23, 31, 12, 1, 2, 3
	 *
	 * @param xyz       单元坐标数组
	 * @param mat       5x5 本构矩阵，下标排列顺序为11, 22, 33, 23, 31, 12, 1, 2, 3
	 * @param tp        问题类型
	 * @param K 得到的单元刚度矩阵
	 */
	public static void elementPEStiff(int[] ele, double[][] xyz, double[][] mat, BVPType tp, double[][] K) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
//		System.out.println(Arrays.toString(x));
//		System.out.println(Arrays.toString(y));
//		System.out.println(Arrays.toString(z));
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] dz = new double[8];
		double[] v = new double[8];
		double xi, et, zt, detJ, wxarea;
		double mE, mNu;
		double[][] KK = new double[32][32];
		double[][] B = new double[9][32];
		MVO.fill(K, 0.0);
//		System.out.println(MVO.toString(mat));
		switch (tp) {
			case CONSTITUTE_MATRIX_COMMON:
				for (int i = 0; i < 2; i++) {
					xi = GP2[i];
					for (int j = 0; j < 2; j++) {
						et = GP2[j];
						for (int k = 0; k < 2; k++) {
							zt = GP2[k];
							detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
//							System.out.println(Arrays.toString(dx));
//							System.out.println(Arrays.toString(dy));
//							System.out.println(Arrays.toString(dz));
							wxarea = GW2[i] * GW2[j] * GW2[k] * detJ;
							for (int ii = 0; ii < 8; ii++) {// 11, 22, 33, 23, 31, 12
								B[0][4 * ii + 0] = dx[ii];
								B[1][4 * ii + 0] = dy[ii];
								B[2][4 * ii + 0] = dz[ii];
								B[3][4 * ii + 1] = dx[ii];
								B[4][4 * ii + 2] = dy[ii];
								B[5][4 * ii + 3] = dz[ii];
//								B[6][4 * ii + 1] = dy[ii]; B[6][4 * ii + 2] = dx[ii];
//								B[7][4 * ii + 2] = dz[ii]; B[7][4 * ii + 3] = dy[ii];
//								B[8][4 * ii + 1] = dz[ii]; B[8][4 * ii + 3] = dx[ii];
								B[6][4 * ii + 2] = dz[ii]; B[6][4 * ii + 3] = dy[ii];
								B[7][4 * ii + 1] = dz[ii]; B[7][4 * ii + 3] = dx[ii];
								B[8][4 * ii + 1] = dy[ii]; B[8][4 * ii + 2] = dx[ii];
							}
//							System.out.println("B: ");
//							System.out.println(MVO.toString(B));
//							System.out.println("wxarea:" + wxarea);
							mul(matmul(transpose(B), matmul(mat, B)), wxarea, KK);
//							System.out.println("KK: ");
//							System.out.println(MVO.toString(KK));
//							System.exit(1);
							for (int ii = 0; ii < 32; ii++) {
								for (int jj = 0; jj < 32; jj++) {
									K[ii][jj] += KK[ii][jj];
								}
							}
						}
					}
				}
				break;
			default:
				throw new IllegalArgumentException("Wrong BVPType!");
		}
//		System.out.println(MVO.toString(K));
//		System.exit(1);
	}

	public void assmebleElectricalStiff(double[][][] constMatCoef, BVPType tp, SRMatrix A, int... domainLabel) {
		int[] range = {0};
		double[][] eleMatrix = new double[8][8];
		EleMVecFunc elementStiff = FES3H81::elementStiff;
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			int eleLabel = mesh.getElementLabel(i);
			for(int k = 0, n = domainLabel.length; k < n; k++) {
				if(eleLabel == domainLabel[k]) {
					mesh.getCoordPerElement(ele, coordRegion);
					elementStiff.action(ele, coordRegion, constMatCoef[0], tp, eleMatrix);
					int[] dof = getDofFromDofs(ele, range);
					A.assemble(dof, eleMatrix);
				}
			}
		}
	}

	public void assmebleElasticStiff(double[][][] constMatCoef, BVPType tp, SRMatrix A, int... domainLabel) {
		int[] range = {1, 2, 3};
		double[][] eleMatrix = new double[24][24];
		EleMVecFunc elementStiff = FES3H82::elementElasticityStiff;
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			int eleLabel = mesh.getElementLabel(i);
			for(int k = 0, n = domainLabel.length; k < n; k++) {
				if(eleLabel == domainLabel[k]) {
					mesh.getCoordPerElement(ele, coordRegion);
					elementStiff.action(ele, coordRegion, constMatCoef[0], tp, eleMatrix);
					int[] dof = getDofFromDofs(ele, range);
					A.assemble(dof, eleMatrixRegion);
				}
			}
		}
	}

	public static void elementPEMass(int[] ele, double[][] xyz, double[][] r, BVPType tp, double[][] K) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] dz = new double[8];
		double[] v = new double[8];
		double xi, et, zt, detJ, wxarea, c1, c2;
		MVO.fill(K, 0.0);
		switch (tp) {
			case CONSIST_MASS:
				for (int i = 0; i < 2; i++) {
					xi = GP2[i];
					for (int j = 0; j < 2; j++) {
						et = GP2[j];
						for (int k = 0; k < 2; k++) {
							zt = GP2[k];
							detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
							c1 = MVO.dot_product(r[0], v);
							c2 = MVO.dot_product(r[1], v);
							wxarea = GW2[i] * GW2[j] * GW2[k] * detJ;
							for (int si = 0; si < 8; si++) {
								for (int sj = 0; sj < 8; sj++) {
									K[4 * si + 0][4 * sj + 0] += v[si] * v[sj] * wxarea * c1;
									K[4 * si + 1][4 * sj + 1] += v[si] * v[sj] * wxarea * c2;
									K[4 * si + 2][4 * sj + 2] += v[si] * v[sj] * wxarea * c2;
									K[4 * si + 3][4 * sj + 3] += v[si] * v[sj] * wxarea * c2;
								}
							}
						}
					}
				}
				break;
			case LUMP_MASS:
				for (int i = 0; i < 2; i++) {
					xi = GP2[i];
					for (int j = 0; j < 2; j++) {
						et = GP2[j];
						for (int k = 0; k < 2; k++) {
							zt = GP2[k];
							detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
							c1 = MVO.dot_product(r[0], v);
							c2 = MVO.dot_product(r[0], v);
							wxarea = GW2[i] * GW2[j] * GW2[k] * detJ;
							for (int si = 0; si < 8; si++) {
								K[4 * si + 0][4 * si + 0] += v[si] * wxarea * c1;
								K[4 * si + 1][3 * si + 1] += v[si] * wxarea * c2;
								K[4 * si + 2][3 * si + 2] += v[si] * wxarea * c2;
								K[4 * si + 3][3 * si + 3] += v[si] * wxarea * c2;
							}
						}
					}
				}
				break;
			default:
				throw new IllegalArgumentException("Illegal BVPType!");
		}
//		System.out.println(MVO.toString(K));
//		System.exit(1);
	}

	public static void elementPEMassStiff(int[] ele, double[][] xyz, double[][][] r, BVPType[] tp, double[][][] K) {
		double[] x = MVO.col(xyz, 0), y = MVO.col(xyz, 1), z = MVO.col(xyz, 2);
		double[] v = new double[8], dx = new double[8], dy = new double[8], dz = new double[8];
		double xi, et, zt, detJ, wxarea;
		double c1, c2;
		MVO.fill(K[0], 0.0);
		MVO.fill(K[1], 0.0);
		double[][] KK = new double[32][32];
		double[][] B = new double[9][32];
		double[][] ConstituteMatrix = r[1];
		for (int i = 0; i < 2; i++) {
			xi = GP2[i];
			for (int j = 0; j < 2; j++) {
				et = GP2[j];
				for (int k = 0; k < 2; k++) {
					zt = GP2[k];
					detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
					c1 = MVO.dot_product(r[0][0], v);
					c2 = MVO.dot_product(r[0][1], v);
					wxarea = GW2[i] * GW2[j] * GW2[k] * detJ;
					for (int ii = 0; ii < 8; ii++) {// 11, 22, 33, 23, 31, 12
						B[0][4 * ii + 0] = dx[ii];
						B[1][4 * ii + 0] = dy[ii];
						B[2][4 * ii + 0] = dz[ii];
						B[3][4 * ii + 1] = dx[ii];
						B[4][4 * ii + 2] = dy[ii];
						B[5][4 * ii + 3] = dz[ii];
//						B[6][4 * ii + 1] = dy[ii]; B[6][4 * ii + 2] = dx[ii];
//						B[7][4 * ii + 2] = dz[ii]; B[7][4 * ii + 3] = dy[ii];
//						B[8][4 * ii + 1] = dz[ii]; B[8][4 * ii + 3] = dx[ii];
						B[6][4 * ii + 2] = dz[ii]; B[6][4 * ii + 3] = dy[ii];
						B[7][4 * ii + 1] = dz[ii]; B[7][4 * ii + 3] = dx[ii];
						B[8][4 * ii + 1] = dy[ii]; B[8][4 * ii + 2] = dx[ii];
					}
					mul(matmul(transpose(B), matmul(ConstituteMatrix, B)), wxarea, KK);
					for (int si = 0; si < 8; si++) {
						for (int sj = 0; sj < 8; sj++) {
							K[0][4 * si + 0][4 * sj + 0] += v[si] * v[sj] * wxarea * c1;
							K[0][4 * si + 1][4 * sj + 1] += v[si] * v[sj] * wxarea * c2;
							K[0][4 * si + 2][4 * sj + 2] += v[si] * v[sj] * wxarea * c2;
							K[0][4 * si + 3][4 * sj + 3] += v[si] * v[sj] * wxarea * c2;
							K[1][si][sj] += KK[si][sj];
						}
					}
				}
			}
		}
		System.out.println("Mass Matrix:");
		System.out.println(MVO.toString(K[0]));
		System.out.println("Stiff Matrix:");
		System.out.println(MVO.toString(K[1]));
	}
}
